Summary:

Two avenues of analysis:

  1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength^2 and \(\Delta\)frequency^2, accounting for bee size and treatment order
  2. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

Measured variables:

Calculated variables:

Install packages and read data

ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS",'reshape2', 
              'influence.ME', 'sjPlot', "effects", 'visreg', 'viridis')
ipak(packages)
##      ggplot2          car         lme4       gsheet         MASS 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
##     reshape2 influence.ME       sjPlot      effects       visreg 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
##      viridis 
##         TRUE
# read in data -- google sheet called "Final Bee Resp Data"
url = 'https://docs.google.com/spreadsheets/d/1wT-QxSJJElhhJcIXlg2hDpFKHbLyFGuqNNj2iYvA8Vo/edit?usp=sharing'
bdta <- data.frame(gsheet2tbl(url))
summary(bdta)
##     BeeID               order      Treatment            Mstarved     
##  Length:60          Min.   :1.0   Length:60          Min.   :0.0767  
##  Class :character   1st Qu.:1.0   Class :character   1st Qu.:0.1213  
##  Mode  :character   Median :1.5   Mode  :character   Median :0.1378  
##                     Mean   :1.5                      Mean   :0.1411  
##                     3rd Qu.:2.0                      3rd Qu.:0.1670  
##                     Max.   :2.0                      Max.   :0.1909  
##        M2               MF             Itspan               S            
##  Min.   :0.1062   Min.   :0.0969   Min.   :0.003860   Min.   :4.010e-05  
##  1st Qu.:0.1808   1st Qu.:0.1762   1st Qu.:0.004400   1st Qu.:5.590e-05  
##  Median :0.2191   Median :0.2107   Median :0.004620   Median :6.130e-05  
##  Mean   :0.2245   Mean   :0.2154   Mean   :0.004594   Mean   :6.141e-05  
##  3rd Qu.:0.2646   3rd Qu.:0.2576   3rd Qu.:0.004870   3rd Qu.:6.750e-05  
##  Max.   :0.3550   Max.   :0.3440   Max.   :0.005180   Max.   :7.810e-05  
##       MetR             freq            amp               L           
##  Min.   : 5.507   Min.   :163.9   Min.   : 96.88   Min.   :0.008691  
##  1st Qu.: 8.935   1st Qu.:173.0   1st Qu.:114.16   1st Qu.:0.010570  
##  Median :10.666   Median :179.1   Median :125.96   Median :0.011017  
##  Mean   :10.600   Mean   :178.9   Mean   :123.31   Mean   :0.010911  
##  3rd Qu.:12.297   3rd Qu.:185.3   3rd Qu.:133.61   3rd Qu.:0.011444  
##  Max.   :16.475   Max.   :194.7   Max.   :144.78   Max.   :0.012118

Calculate variables

bdta <- within(bdta, {
     MT = (M2 + MF)/2
     load = MT - Mstarved
     perLoad = load / Mstarved * 100
     arcL = (0.75 * L) * (amp * (pi / 180))
     U = arcL * freq * 2
     freq2 = freq^2
     arcL2 = arcL^2
     frce = U^2 * S
     
     # convert to factor variables
     order = as.factor(as.character(order))
     Treatment = as.factor(as.character(Treatment))
     BeeID = as.factor(as.character(BeeID))
})

Create a new dataframe that calculates the changes for each bee

# convert from long to wide
newDF <- data.frame()
colsTocalc = c("order", "M2", "MF", "MetR", "freq", "amp", "load", "MT", "perLoad", "frce", "arcL2", "freq2", "U", "arcL")
for(varb in colsTocalc){
     data_wide <- dcast(bdta, BeeID + Mstarved + S + Itspan + L ~ 
                             Treatment, value.var=c(varb))
     colnames(data_wide)[6:7] = paste(varb, colnames(data_wide)[6:7], sep = "_")
     if(varb == colsTocalc[1]){
       newDF <- data_wide   
     }
     else newDF <- merge(newDF, data_wide, all.y = TRUE)    
}

head(newDF)
##   BeeID Mstarved        S  Itspan          L order_H order_L   M2_H   M2_L
## 1   E01   0.1577 6.46e-05 0.00464 0.01116843       2       1 0.2634 0.1895
## 2   E03   0.1732 6.92e-05 0.00487 0.01144364       1       2 0.3106 0.2105
## 3   E05   0.1755 7.06e-05 0.00497 0.01168730       2       1 0.2776 0.2006
## 4   E09   0.1135 4.88e-05 0.00440 0.00986884       1       2 0.2282 0.1425
## 5   E10   0.1269 5.91e-05 0.00462 0.01060036       1       2 0.2166 0.1460
## 6   E11   0.1350 5.88e-05 0.00464 0.01070018       1       2 0.2223 0.1599
##     MF_H   MF_L   MetR_H    MetR_L   freq_H   freq_L    amp_H     amp_L
## 1 0.2559 0.1812 12.54241 10.393483 186.9973 179.3308 128.7145 111.49875
## 2 0.3079 0.2023 11.78387  9.214622 174.9374 167.8412 144.7784 119.98953
## 3 0.2742 0.1975 12.20102  8.981412 178.8896 166.2385 125.6778 116.22053
## 4 0.2084 0.1390  9.47114  7.581104 191.8348 190.3535 139.3246 103.30703
## 5 0.2099 0.1422 11.03371  7.159840 186.1762 169.5451 128.8257  96.88003
## 6 0.2179 0.1579 10.23580  6.806696 183.0259 168.2223 128.7916 116.18262
##    load_H  load_L    MT_H    MT_L perLoad_H perLoad_L      frce_H
## 1 0.10195 0.02765 0.25965 0.18535  64.64807  17.53329 0.003199482
## 2 0.13605 0.03320 0.30925 0.20640  78.55081  19.16859 0.003984230
## 3 0.10040 0.02355 0.27590 0.19905  57.20798  13.41880 0.003340857
## 4 0.10480 0.02725 0.21830 0.14075  92.33480  24.00881 0.002327019
## 5 0.08635 0.01720 0.21325 0.14410  68.04571  13.55398 0.002618299
## 6 0.08510 0.02390 0.22010 0.15890  63.03704  17.70370 0.002563877
##        frce_L      arcL2_H      arcL2_L  freq2_H  freq2_L      U_H
## 1 0.002208025 0.0003540922 0.0002657061 34967.99 32159.55 7.037583
## 2 0.002519159 0.0004703414 0.0003230668 30603.08 28170.67 7.587858
## 3 0.002467172 0.0003696773 0.0003161342 32001.50 27635.24 6.879020
## 4 0.001259710 0.0003239404 0.0001781022 36800.61 36234.44 6.905419
## 5 0.001228018 0.0003195386 0.0001807120 34661.58 28745.53 6.656039
## 6 0.001762569 0.0003254128 0.0002648146 33498.49 28298.73 6.603283
##        U_L     arcL_H     arcL_L
## 1 5.846363 0.01881734 0.01630050
## 2 6.033576 0.02168736 0.01797406
## 3 5.911496 0.01922699 0.01778016
## 4 5.080722 0.01799834 0.01334549
## 5 4.558361 0.01787564 0.01344292
## 6 5.475004 0.01803920 0.01627313

Calculate \(\Delta\) variables

newDF <- within(newDF, {
     deltaPercLoad = perLoad_H - perLoad_L
     avgPercLoad = (perLoad_H + perLoad_L) / 2
     deltaMetR = MetR_H - MetR_L
     deltaFrq2 = freq2_H - freq2_L
     deltaArcL = arcL_H - arcL_L
     deltaArcL2 = arcL2_H - arcL2_L
     deltaFreq2Perc = deltaFrq2 / deltaPercLoad
     deltaLoad = scale(load_H - load_L, center = TRUE, scale = FALSE)
     dLoad_nonCent = load_H - load_L
     deltaLoad2 <- deltaLoad^2
})


plot(newDF$deltaArcL2 ~ newDF$deltaFrq2)

# compare with Susie's calculations



b2 <- read.csv("~/Desktop/b2.csv")
b2 <- b2[!is.na(b2$deltaperload), ]
b2 <- b2[order(b2$BeeID, b2$Treatment), ]

plot(b2$AverageperLoad, newDF$avgPercLoad)
abline(0,1)

plot(b2$deltaload, newDF$deltaLoad)
abline(0,1)

Use PCA to combine bee sizes into a single predictor

# principle components
aa = prcomp(newDF[, c("Mstarved", "Itspan", "S", "L")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
##                           PC1     PC2     PC3     PC4
## Standard deviation     1.9422 0.37234 0.22746 0.19332
## Proportion of Variance 0.9431 0.03466 0.01293 0.00934
## Cumulative Proportion  0.9431 0.97772 0.99066 1.00000
biplot(aa) # shows that all three size measurement are correlated

# note, I changed the signs of the predictions so that higher PC1 values 
# correspond to bigger bees
p1 = -predict(aa)[,1] 

# add PC1 scores to dataset
newDF$size_pc1 = p1
newDF$size_pc1_2 = newDF$size_pc1^2

# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(newDF[, c("Mstarved", "Itspan", "S",  "L",  "size_pc1")])

1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength accounting for bee size and treatment order

This is the model selection procedure that was used: 1. Fit a large model with all two-way interactions and squared terms 2. Remove non-significant predictors, starting with interactions, squared terms, and then main effects

# reformat order so that it is more interpretable
library(plyr)
newDF$order_1 <- mapvalues(newDF$order_H, from = c(2, 1), to = c("loadedSecond", "loadedFirst"))

# fit full model
m1 <- lm(deltaArcL2 ~  (deltaLoad +  size_pc1 + order_1)^2  + 
              size_pc1_2  + deltaLoad2, data = newDF)
summary(m1)
## 
## Call:
## lm(formula = deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + 
##     size_pc1_2 + deltaLoad2, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.326e-05 -1.951e-05  3.244e-06  1.635e-05  4.401e-05 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.818e-05  1.300e-05   7.551 2.05e-07 ***
## deltaLoad                      1.859e-03  1.178e-03   1.578    0.130    
## size_pc1                      -1.804e-05  1.263e-05  -1.428    0.168    
## order_1loadedSecond            7.310e-07  1.978e-05   0.037    0.971    
## size_pc1_2                    -3.660e-06  4.905e-06  -0.746    0.464    
## deltaLoad2                     3.807e-02  4.366e-02   0.872    0.393    
## deltaLoad:size_pc1            -2.300e-06  8.387e-04  -0.003    0.998    
## deltaLoad:order_1loadedSecond -6.487e-04  2.225e-03  -0.292    0.774    
## size_pc1:order_1loadedSecond   1.204e-05  2.099e-05   0.574    0.572    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.03e-05 on 21 degrees of freedom
## Multiple R-squared:  0.575,  Adjusted R-squared:  0.413 
## F-statistic: 3.551 on 8 and 21 DF,  p-value: 0.009387
m2 <- update(m1, .~. - deltaLoad:size_pc1)
anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 + 
##     deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1 + size_pc1:order_1
##   Res.Df        RSS Df   Sum of Sq  F Pr(>F)
## 1     21 1.9281e-08                         
## 2     22 1.9281e-08 -1 -6.9071e-15  0 0.9978
summary(m2)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:order_1 + size_pc1:order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.327e-05 -1.953e-05  3.263e-06  1.634e-05  4.402e-05 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.817e-05  1.252e-05   7.840 8.24e-08 ***
## deltaLoad                      1.860e-03  9.935e-04   1.872   0.0745 .  
## size_pc1                      -1.806e-05  7.886e-06  -2.290   0.0320 *  
## order_1loadedSecond            7.471e-07  1.846e-05   0.040   0.9681    
## size_pc1_2                    -3.672e-06  1.745e-06  -2.104   0.0470 *  
## deltaLoad2                     3.800e-02  3.336e-02   1.139   0.2669    
## deltaLoad:order_1loadedSecond -6.509e-04  2.022e-03  -0.322   0.7505    
## size_pc1:order_1loadedSecond   1.209e-05  1.187e-05   1.019   0.3194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.96e-05 on 22 degrees of freedom
## Multiple R-squared:  0.575,  Adjusted R-squared:  0.4397 
## F-statistic: 4.251 on 7 and 22 DF,  p-value: 0.004168
m3 <- update(m2, .~. - deltaLoad:order_1)

anova(m2, m3)
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1 + size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     size_pc1:order_1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     22 1.9281e-08                             
## 2     23 1.9371e-08 -1 -9.0841e-11 0.1037 0.7505
summary(m3)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + size_pc1:order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.414e-05 -1.856e-05  3.421e-06  1.715e-05  4.272e-05 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   9.768e-05  1.218e-05   8.017 4.13e-08 ***
## deltaLoad                     1.625e-03  6.585e-04   2.467   0.0215 *  
## size_pc1                     -1.728e-05  7.355e-06  -2.349   0.0278 *  
## order_1loadedSecond           3.683e-06  1.573e-05   0.234   0.8170    
## size_pc1_2                   -3.816e-06  1.654e-06  -2.308   0.0304 *  
## deltaLoad2                    4.504e-02  2.467e-02   1.825   0.0809 .  
## size_pc1:order_1loadedSecond  1.003e-05  9.807e-06   1.023   0.3169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.902e-05 on 23 degrees of freedom
## Multiple R-squared:  0.573,  Adjusted R-squared:  0.4615 
## F-statistic: 5.143 on 6 and 23 DF,  p-value: 0.001756
m4 <- update(m3, .~. - size_pc1:order_1)
anova(m3, m4) 
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     23 1.9371e-08                             
## 2     24 2.0253e-08 -1 -8.8175e-10 1.0469 0.3169
summary(m4)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.427e-05 -1.874e-05 -5.000e-09  1.783e-05  4.254e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.024e-04  1.128e-05   9.079 3.13e-09 ***
## deltaLoad            1.639e-03  6.590e-04   2.486   0.0203 *  
## size_pc1            -1.180e-05  5.042e-06  -2.340   0.0279 *  
## order_1loadedSecond -5.411e-07  1.519e-05  -0.036   0.9719    
## size_pc1_2          -2.627e-06  1.178e-06  -2.231   0.0353 *  
## deltaLoad2           2.613e-02  1.636e-02   1.597   0.1233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.905e-05 on 24 degrees of freedom
## Multiple R-squared:  0.5535, Adjusted R-squared:  0.4605 
## F-statistic: 5.951 on 5 and 24 DF,  p-value: 0.001027
m5 <- update(m4, .~. - deltaLoad2)
anova(m4,m5) 
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     24 2.0253e-08                             
## 2     25 2.2406e-08 -1 -2.1529e-09 2.5511 0.1233
summary(m5)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2, 
##     data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.782e-05 -1.748e-05  9.200e-07  1.678e-05  4.278e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.074e-04  1.117e-05   9.623 6.93e-10 ***
## deltaLoad            2.095e-03  6.122e-04   3.422  0.00215 ** 
## size_pc1            -1.388e-05  5.018e-06  -2.767  0.01050 *  
## order_1loadedSecond  2.783e-07  1.565e-05   0.018  0.98595    
## size_pc1_2          -2.037e-06  1.152e-06  -1.768  0.08930 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.994e-05 on 25 degrees of freedom
## Multiple R-squared:  0.5061, Adjusted R-squared:  0.427 
## F-statistic: 6.403 on 4 and 25 DF,  p-value: 0.001086
m6 <- update(m5, .~.  - size_pc1_2)
anova(m6, m5) 
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)  
## 1     26 2.5207e-08                              
## 2     25 2.2406e-08  1 2.8009e-09 3.1251 0.0893 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m6)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.154e-05 -1.781e-05 -4.100e-07  1.849e-05  4.983e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          9.844e-05  1.034e-05   9.525 5.79e-10 ***
## deltaLoad            2.148e-03  6.359e-04   3.378  0.00231 ** 
## size_pc1            -1.321e-05  5.205e-06  -2.539  0.01745 *  
## order_1loadedSecond  3.224e-06  1.618e-05   0.199  0.84364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.114e-05 on 26 degrees of freedom
## Multiple R-squared:  0.4443, Adjusted R-squared:  0.3802 
## F-statistic: 6.929 on 3 and 26 DF,  p-value: 0.001402
m7 <- update(m6, .~. - order_1)
anova(m7, m6)
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1     27 2.5245e-08                            
## 2     26 2.5207e-08  1 3.8477e-11 0.0397 0.8436
summary(m7)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666
anova(m7, update(m7, .~. - deltaLoad)) # p-value for delta load
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ size_pc1
##   Res.Df        RSS Df   Sum of Sq      F    Pr(>F)    
## 1     27 2.5245e-08                                    
## 2     28 4.5359e-08 -1 -2.0114e-08 21.512 8.055e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m7, update(m7, .~. - size_pc1)) # p-value for size
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad
##   Res.Df        RSS Df   Sum of Sq      F   Pr(>F)   
## 1     27 2.5245e-08                                  
## 2     28 3.4619e-08 -1 -9.3737e-09 10.025 0.003808 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m7)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666
# refit model with non-centered version of deltaLoad
# this model will have a different intercept, but same p-values for slopes
summary(update(m7, .~. - deltaLoad + dLoad_nonCent)) # final model for paper
## 
## Call:
## lm(formula = deltaArcL2 ~ size_pc1 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -5.822e-05  3.460e-05  -1.683  0.10397    
## size_pc1      -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## dLoad_nonCent  2.059e-03  4.439e-04   4.638 8.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666

model diagnostics

par(mfrow = c(2,2))
plot(m7, which = 1:4) # no glaring violations, though there are a few fairly influential observations

par(mfrow = c(1,1))

car::vif(m7) 
## deltaLoad  size_pc1 
##  1.841066  1.841066

Model visualization

summary(m7)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666
# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(m7, type = 'partial')[, "deltaLoad"]

# plot partial residuals with base R plotting
plot(x = newDF$deltaLoad, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) 
## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.584e-04  2.531e-05  -6.257 9.20e-07 ***
## newDF$dLoad_nonCent  2.059e-03  3.213e-04   6.409 6.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.003e-05 on 28 degrees of freedom
## Multiple R-squared:  0.5946, Adjusted R-squared:  0.5801 
## F-statistic: 41.07 on 1 and 28 DF,  p-value: 6.137e-07
# this is what the raw data look like
plot(x = newDF$deltaLoad, y = newDF$deltaFrq2)

# this package plots the partial residuals
crPlot(m7, variable = "deltaLoad")

# plot with ggplot2
theme_set(theme_classic()) 

# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaArcL, color = deltaLoad)) + 
     geom_point()

ggplot(newDF, aes(x= deltaLoad, y = deltaArcL, color = size_pc1)) + 
     geom_point()

# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) + 
     geom_point() + 
     labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on arcL^2 \n while subtracing affect of bee size") + 
     stat_smooth(method = 'lm', se = FALSE)

partialResidDeltaLoad <- data.frame(deltaLoad = newDF$dLoad_nonCent, partResArcL2 = y )

Takeaways from model 1:

  1. A higher deltaload (i.e. a relatively larger “high” load) causes bees to have a larger change in arclength squared. Another way to say this is that increasing the change in load by 1 gram causes an increase in the change in arclength^2 by 0.002059 degrees.
  2. The larger the bee, the lower the change in arcLength^2 (association, rather than causation), while holding the change in load constant.
  3. We do not have enough evidence to notice any non-linear relationships in this model.

See how \(\Delta\)load affects \(\Delta\)frequency, accounting for bee size and treatment order

# fit full model
m1 <- lm(deltaFrq2 ~  (deltaLoad +  size_pc1 + order_1)^2  + 
              size_pc1_2  + deltaLoad2, data = newDF)
summary(m1)
## 
## Call:
## lm(formula = deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + 
##     size_pc1_2 + deltaLoad2, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3981.4 -1255.9  -170.3  1165.6  3552.9 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4195.97     909.39   4.614  0.00015 ***
## deltaLoad                     -106368.80   82390.69  -1.291  0.21073    
## size_pc1                          666.06     883.26   0.754  0.45917    
## order_1loadedSecond             -4376.19    1383.48  -3.163  0.00469 ** 
## size_pc1_2                       -120.74     343.06  -0.352  0.72839    
## deltaLoad2                    -469149.49 3053199.27  -0.154  0.87935    
## deltaLoad:size_pc1               6902.45   58657.40   0.118  0.90744    
## deltaLoad:order_1loadedSecond   62840.01  155625.35   0.404  0.69045    
## size_pc1:order_1loadedSecond       54.39    1467.95   0.037  0.97079    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2119 on 21 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.3162 
## F-statistic: 2.676 on 8 and 21 DF,  p-value: 0.03372
car::vif(m1) # some serious multicollinearity
##          deltaLoad           size_pc1            order_1 
##          13.203325          19.003708           3.182227 
##         size_pc1_2         deltaLoad2 deltaLoad:size_pc1 
##          18.740891          10.117446          25.115087 
##  deltaLoad:order_1   size_pc1:order_1 
##          11.384208          23.301468
m2 <- update(m1, .~. - size_pc1:order_1)
anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 + 
##     deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:size_pc1 + deltaLoad:order_1
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     21 94310333                           
## 2     22 94316499 -1   -6165.2 0.0014 0.9708
summary(m2)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:size_pc1 + deltaLoad:order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3973.1 -1259.6  -176.2  1149.9  3556.7 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4205.4      852.7   4.932 6.21e-05 ***
## deltaLoad                     -108212.1    64166.1  -1.686  0.10584    
## size_pc1                          695.7      367.7   1.892  0.07176 .  
## order_1loadedSecond             -4386.5     1324.0  -3.313  0.00316 ** 
## size_pc1_2                       -109.9      176.0  -0.625  0.53864    
## deltaLoad2                    -417980.4  2660489.9  -0.157  0.87659    
## deltaLoad:size_pc1               5130.1    33170.6   0.155  0.87850    
## deltaLoad:order_1loadedSecond   66238.1   122843.3   0.539  0.59516    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2071 on 22 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.3472 
## F-statistic: 3.203 on 7 and 22 DF,  p-value: 0.01701
m3 <- update(m2, .~. - deltaLoad:size_pc1)

anova(m2, m3)
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:size_pc1 + deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     22 94316499                           
## 2     23 94419044 -1   -102545 0.0239 0.8785
summary(m3)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3910.3 -1244.0  -168.5  1173.4  3503.8 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4190.99     829.34   5.053 4.09e-05 ***
## deltaLoad                     -108478.92   62767.13  -1.728  0.09734 .  
## size_pc1                          702.89     356.93   1.969  0.06108 .  
## order_1loadedSecond             -4434.40    1259.70  -3.520  0.00184 ** 
## size_pc1_2                        -88.67     107.58  -0.824  0.41829    
## deltaLoad2                    -215246.26 2265485.45  -0.095  0.92513    
## deltaLoad:order_1loadedSecond   61643.10  116639.62   0.528  0.60222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2026 on 23 degrees of freedom
## Multiple R-squared:  0.5042, Adjusted R-squared:  0.3749 
## F-statistic: 3.899 on 6 and 23 DF,  p-value: 0.00785
m4 <- update(m3, .~. - deltaLoad:order_1)
anova(m3, m4) 
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     23 94419044                           
## 2     24 95565635 -1  -1146591 0.2793 0.6022
summary(m4)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3667.1 -1308.7   -52.9  1317.9  3482.6 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.330e+03  7.749e+02   5.587 9.48e-06 ***
## deltaLoad           -8.589e+04  4.527e+04  -1.897 0.069889 .  
## size_pc1             7.352e+02  3.463e+02   2.123 0.044282 *  
## order_1loadedSecond -4.794e+03  1.044e+03  -4.594 0.000117 ***
## size_pc1_2          -5.195e+01  8.090e+01  -0.642 0.526837    
## deltaLoad2          -1.250e+06  1.124e+06  -1.112 0.277131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1995 on 24 degrees of freedom
## Multiple R-squared:  0.4982, Adjusted R-squared:  0.3937 
## F-statistic: 4.766 on 5 and 24 DF,  p-value: 0.00364
m5 <- update(m4, .~. - size_pc1_2)
anova(m4,m5) 
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     24 95565635                           
## 2     25 97207799 -1  -1642165 0.4124 0.5268
summary(m5)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3591.2 -1388.0  -154.6  1391.9  3336.4 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             4166.0      723.2   5.760 5.29e-06 ***
## deltaLoad             -80714.2    44020.7  -1.834 0.078654 .  
## size_pc1                 732.5      342.2   2.141 0.042251 *  
## order_1loadedSecond    -4719.5     1024.9  -4.605 0.000104 ***
## deltaLoad2          -1475818.2  1054446.2  -1.400 0.173916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1972 on 25 degrees of freedom
## Multiple R-squared:  0.4896, Adjusted R-squared:  0.4079 
## F-statistic: 5.995 on 4 and 25 DF,  p-value: 0.00159
m6 <- update(m5, .~.  - deltaLoad2)
anova(m6, m5) 
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1     26 104824695                           
## 2     25  97207799  1   7616896 1.9589 0.1739
summary(m6)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3735.3      666.5   5.605 6.87e-06 ***
## deltaLoad           -105594.7    41007.6  -2.575 0.016065 *  
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244
m7 <- update(m6, .~. - size_pc1)
anova(m7, m6) # p-values for size
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1     27 131380648                              
## 2     26 104824695  1  26555953 6.5868 0.01638 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m8 <- update(m6, .~. - deltaLoad)
anova(m6, m8) # p-value for deltaLoad
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ size_pc1 + order_1
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1     26 104824695                              
## 2     27 131557512 -1 -26732816 6.6306 0.01607 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m9 <- update(m6, .~. - order_1)
anova(m6, m9) # p-value for order
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1     26 104824695                                  
## 2     27 187214161 -1 -82389465 20.435 0.0001192 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit final model with non-centered deltaLoad
m10 <- update(m6, .~. - deltaLoad + dLoad_nonCent)
summary(m10) # final model for paper
## 
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           11857.7     3586.6   3.306 0.002766 ** 
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## dLoad_nonCent       -105594.7    41007.6  -2.575 0.016065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

model diagnostics

par(mfrow = c(2,2))
plot(m10, which = 1:4) # no glaring violations

par(mfrow = c(1,1))

car::vif(m10) # vif is a little high
##      size_pc1       order_1 dLoad_nonCent 
##      3.056456      2.017003      3.643394

model visualization

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(m10, type = 'partial')[, "dLoad_nonCent"]

partialResidDeltaLoad <- cbind(partialResidDeltaLoad, y)
colnames(partialResidDeltaLoad)[3] <- "partResFreq2"


write.csv(partialResidDeltaLoad, file = "PartialResidArcL2.csv", row.names = FALSE)

# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) 
## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             8122       1631   4.980 2.93e-05 ***
## newDF$dLoad_nonCent  -105595      20702  -5.101 2.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1935 on 28 degrees of freedom
## Multiple R-squared:  0.4816, Adjusted R-squared:  0.4631 
## F-statistic: 26.02 on 1 and 28 DF,  p-value: 2.106e-05
# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaFrq2)

# this package plots the partial residuals
crPlot(m10, variable = "dLoad_nonCent")

# plot with ggplot2
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaFrq2, color = deltaLoad, shape = order_1)) + 
     geom_point()

ggplot(newDF, aes(x= deltaLoad, y = deltaFrq2, color = size_pc1, shape = order_1)) + 
     geom_point()

# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) + 
     geom_point() + 
     labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on freq^2 \n while subtracing affect of bee size and order") + 
     stat_smooth(method = 'lm', se = FALSE)

# partial residuals for order

y <- residuals(m10, type = 'partial')[, "order_1"]

ggplot(newDF, aes(x= order_1, y = y)) + 
     geom_boxplot() + 
     labs(x = "order", y = "partial residuals for order \n i.e. order effect on freq^2 \n while subtracing affect of bee size and load") + 
     stat_smooth(method = 'lm', se = FALSE)

# holding bee size and delta load constant, if a bee was loaded second, (confusing, huh?) then it would have a much lower freq^2 than if it was loaded first

summary(m10)
## 
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           11857.7     3586.6   3.306 0.002766 ** 
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## dLoad_nonCent       -105594.7    41007.6  -2.575 0.016065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

Takeaways from model 2:

  1. Order, size, and deltaLoad are associated with a change in frequency^2. We find no evidence of non-linear relationships or interactions.
  2. Holding size and order constant, we find that a larger deltaLoad causes a decrease in deltaFrequency^2
  3. Holding other variables constant, an increase in bee size is associated with a larger deltaFrequency^2.
  4. Holding other variables constant, if the bee was loaded second, then they had a lower deltaFrequency^2 than if they were loaded first.
# make some predictions
predDF <- data.frame(size_pc1 = 0, order_1 = c("loadedSecond", "loadedFirst"), dLoad_nonCent = mean(newDF$dLoad_nonCent))

predDF$pred_deltaF2 <- predict(m10, predDF)

predDF
##   size_pc1      order_1 dLoad_nonCent pred_deltaF2
## 1        0 loadedSecond       0.07692     -982.314
## 2        0  loadedFirst       0.07692     3735.337

2. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

mm1 <- lm(deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2, data = newDF)
car::vif(mm1)
##  deltaFrq2 deltaArcL2  deltaLoad   size_pc1 deltaLoad2 
##   1.060059   1.852143   3.844322   2.679591   1.495113
summary(mm1)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + 
##     size_pc1 + deltaLoad2, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11198 -0.64060 -0.01055  0.42457  2.97948 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.941e+00  6.182e-01   3.140  0.00444 ** 
## deltaFrq2    3.543e-04  6.888e-05   5.144 2.89e-05 ***
## deltaArcL2  -2.573e+03  5.900e+03  -0.436  0.66659    
## deltaLoad    2.368e+01  1.937e+01   1.223  0.23337    
## size_pc1    -5.482e-02  1.445e-01  -0.379  0.70776    
## deltaLoad2   4.440e+02  5.114e+02   0.868  0.39381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9233 on 24 degrees of freedom
## Multiple R-squared:  0.5937, Adjusted R-squared:  0.5091 
## F-statistic: 7.014 on 5 and 24 DF,  p-value: 0.0003629
mm2 <- update(mm1, .~. - size_pc1)
anova(mm1, mm2)
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     24 20.460                           
## 2     25 20.582 -1  -0.12268 0.1439 0.7078
summary(mm2)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06884 -0.60009 -0.03648  0.36122  3.09536 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.821e+00  5.223e-01   3.487  0.00182 ** 
## deltaFrq2    3.571e-04  6.730e-05   5.306 1.69e-05 ***
## deltaArcL2  -1.543e+03  5.147e+03  -0.300  0.76682    
## deltaLoad    1.788e+01  1.168e+01   1.530  0.13849    
## deltaLoad2   4.891e+02  4.888e+02   1.001  0.32665    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9074 on 25 degrees of freedom
## Multiple R-squared:  0.5913, Adjusted R-squared:  0.5259 
## F-statistic: 9.041 on 4 and 25 DF,  p-value: 0.0001166
mm3 <- update(mm2, .~. - deltaArcL2)
anova(mm3, mm2)
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 20.656                           
## 2     25 20.582  1   0.07399 0.0899 0.7668
summary(mm3)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2, 
##     data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05850 -0.58135 -0.00991  0.37733  3.10956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.682e+00  2.327e-01   7.228 1.12e-07 ***
## deltaFrq2   3.560e-04  6.602e-05   5.393 1.19e-05 ***
## deltaLoad   1.667e+01  1.077e+01   1.548    0.134    
## deltaLoad2  4.422e+02  4.550e+02   0.972    0.340    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8913 on 26 degrees of freedom
## Multiple R-squared:  0.5898, Adjusted R-squared:  0.5425 
## F-statistic: 12.46 on 3 and 26 DF,  p-value: 3.065e-05
mm4 <- update(mm3, .~. - deltaLoad2)
anova(mm3, mm4)
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 20.656                           
## 2     27 21.407 -1  -0.75053 0.9447   0.34
summary(mm4)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.824e+00  1.808e-01  10.086 1.18e-10 ***
## deltaFrq2   3.451e-04  6.498e-05   5.310 1.32e-05 ***
## deltaLoad   2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06
mm5 <- update(mm4, .~. - deltaLoad) 
anova(mm5, mm4) # p-value for load
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     28 25.352                              
## 2     27 21.407  1    3.9448 4.9754 0.03421 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mm6 <- update(mm4, .~. - deltaFrq2)
anova(mm4, mm6) # p-value for deltafrq2
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad
## Model 2: deltaMetR ~ deltaLoad
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     27 21.407                                  
## 2     28 43.766 -1   -22.359 28.201 1.324e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit model with non-centered load
mm7 <- update(mm4, .~. - deltaLoad + dLoad_nonCent)

summary(mm7) # final model for paper
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

model diagnostics

car::vif(mm7)
##     deltaFrq2 dLoad_nonCent 
##      1.014369      1.014369
par(mfrow = c(2,2))
plot(mm7, which = 1:4) # no glaring violations

par(mfrow = c(1,1))

model visualization

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(mm7, type = 'partial')[, "dLoad_nonCent"]

# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) 
## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          -1.6463     0.7371  -2.233   0.0337 *
## newDF$dLoad_nonCent  21.4030     9.3554   2.288   0.0299 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8744 on 28 degrees of freedom
## Multiple R-squared:  0.1575, Adjusted R-squared:  0.1274 
## F-statistic: 5.234 on 1 and 28 DF,  p-value: 0.02991
summary(mm7)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06
# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaMetR)

# this package plots the partial residuals
crPlot(mm7, variable = "dLoad_nonCent")

crPlot(mm7, variable = "deltaFrq2")

summary(mm7) # final model for paper
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

Takeaways for model 3:

  1. Deltafreq^2 and deltaLoad are both associated with deltaMetR.
  2. Holding deltafreq constant (not necessarily at 0), an increase in deltaload is associated with an increase in deltametabolicRate
  3. Holding deltaLoad constant (not holding load constant, and again, not holding deltaload necessarily at 0), an increase in deltaFreq^2 is associated with an increase in deltaMetabolic rate.
  4. We found no evidence to suggest that a change in deltaArcL^2 is associated with deltametabolicRate. (this is not saying that arcLength doesn’t affect metabolic rate).

Refref: Do regression for Average percent loading vs. delta metabolic rate / 1% load etc. It’s the last page of the document Susie sent (3 different response variables).

Covariates: avg % load, order, bee size.

(DeltaMetRate / deltaPerLoad) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {dmr_dpl = deltaMetR / deltaPercLoad})

mod1 <- lm(dmr_dpl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1)
## 
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.025580 -0.011566 -0.004357  0.009067  0.043391 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.1078453  0.0158826   6.790 3.31e-07 ***
## avgPercLoad         -0.0010569  0.0002784  -3.796 0.000794 ***
## order_1loadedSecond -0.0131290  0.0066538  -1.973 0.059196 .  
## size_pc1             0.0014674  0.0017912   0.819 0.420121    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01773 on 26 degrees of freedom
## Multiple R-squared:  0.4844, Adjusted R-squared:  0.4249 
## F-statistic: 8.141 on 3 and 26 DF,  p-value: 0.0005502
mod2 <- update(mod1, .~. - size_pc1)
anova(mod1, mod2) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: dmr_dpl ~ avgPercLoad + order_1 + size_pc1
## Model 2: dmr_dpl ~ avgPercLoad + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     26 0.0081723                             
## 2     27 0.0083833 -1 -0.00021093 0.6711 0.4201
mod3 <- update(mod2, .~. - order_1)
anova(mod2, mod3) # remove order
## Analysis of Variance Table
## 
## Model 1: dmr_dpl ~ avgPercLoad + order_1
## Model 2: dmr_dpl ~ avgPercLoad
##   Res.Df       RSS Df  Sum of Sq      F  Pr(>F)  
## 1     27 0.0083833                               
## 2     28 0.0094654 -1 -0.0010822 3.4854 0.07281 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod3) # final model for paper
## 
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.026590 -0.011771 -0.000679  0.011114  0.046931 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.1082986  0.0158566   6.830 2.02e-07 ***
## avgPercLoad -0.0011884  0.0002735  -4.345 0.000165 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01839 on 28 degrees of freedom
## Multiple R-squared:  0.4028, Adjusted R-squared:  0.3814 
## F-statistic: 18.88 on 1 and 28 DF,  p-value: 0.0001654
# get p-values for avgPercLoad and order
anova(mod2, update(mod3, .~. - avgPercLoad)) # p-value for avg perc load
## Analysis of Variance Table
## 
## Model 1: dmr_dpl ~ avgPercLoad + order_1
## Model 2: dmr_dpl ~ 1
##   Res.Df       RSS Df  Sum of Sq      F    Pr(>F)    
## 1     27 0.0083833                                   
## 2     29 0.0158487 -2 -0.0074654 12.022 0.0001846 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plot(I(deltaMetR)  ~ avgPercLoad, data = newDF)
# plot(I(deltaMetR / avgPercLoad)  ~ deltaMetR, data = newDF)

car::scatterplotMatrix(newDF[, c("avgPercLoad", "deltaMetR", "size_pc1")])

car::vif(mod1)  # looks fine
## avgPercLoad     order_1    size_pc1 
##    1.114454    1.051699    1.116654
par(mfrow = c(2,2))
plot(mod3) # possible nonlinear trend in residuals

plot(mod3, which = 4)

# update model to add a non-linear term
newDF <- within(newDF, {avgPercLoad_cent = as.numeric(scale(newDF$avgPercLoad, center = TRUE, scale = FALSE))})
newDF$apl_cent2 <- with(newDF, avgPercLoad_cent^2)

m11 <- lm(dmr_dpl  ~ avgPercLoad_cent  + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m11)
## 
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + 
##     size_pc1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.026372 -0.011430 -0.004154  0.009170  0.044752 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.722e-02  5.369e-03   8.795    4e-09 ***
## avgPercLoad_cent    -1.070e-03  2.860e-04  -3.741  0.00096 ***
## apl_cent2            7.051e-06  2.112e-05   0.334  0.74123    
## order_1loadedSecond -1.373e-02  7.006e-03  -1.960  0.06125 .  
## size_pc1             1.494e-03  1.824e-03   0.819  0.42065    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01804 on 25 degrees of freedom
## Multiple R-squared:  0.4866, Adjusted R-squared:  0.4045 
## F-statistic: 5.925 on 4 and 25 DF,  p-value: 0.0017
car::vif(m11) # vif is much better with centered variables
## avgPercLoad_cent        apl_cent2          order_1         size_pc1 
##         1.135857         1.108436         1.126036         1.118753
m11a <- update(m11, .~. - size_pc1)
anova(m11, m11a) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     25 0.0081360                             
## 2     26 0.0083542 -1 -0.00021818 0.6704 0.4206
m11b <- update(m11a, .~. - apl_cent2)
anova(m11a, m11b) # drop squared term
## Analysis of Variance Table
## 
## Model 1: dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1
## Model 2: dmr_dpl ~ avgPercLoad_cent + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     26 0.0083542                             
## 2     27 0.0083833 -1 -2.9042e-05 0.0904 0.7661
summary(m11b)
## 
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad_cent + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030196 -0.012515 -0.003418  0.009765  0.040799 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.0474326  0.0047309  10.026 1.34e-10 ***
## avgPercLoad_cent    -0.0011242  0.0002644  -4.253 0.000226 ***
## order_1loadedSecond -0.0121420  0.0065038  -1.867 0.072809 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01762 on 27 degrees of freedom
## Multiple R-squared:  0.471,  Adjusted R-squared:  0.4319 
## F-statistic: 12.02 on 2 and 27 DF,  p-value: 0.0001846
m11c <- update(m11b, .~. -order_1 )
anova(m11b, m11c) # remove order_1
## Analysis of Variance Table
## 
## Model 1: dmr_dpl ~ avgPercLoad_cent + order_1
## Model 2: dmr_dpl ~ avgPercLoad_cent
##   Res.Df       RSS Df  Sum of Sq      F  Pr(>F)  
## 1     27 0.0083833                               
## 2     28 0.0094654 -1 -0.0010822 3.4854 0.07281 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m11c)
## 
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad_cent, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.026590 -0.011771 -0.000679  0.011114  0.046931 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.0409569  0.0033568  12.201 1.01e-12 ***
## avgPercLoad_cent -0.0011884  0.0002735  -4.345 0.000165 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01839 on 28 degrees of freedom
## Multiple R-squared:  0.4028, Adjusted R-squared:  0.3814 
## F-statistic: 18.88 on 1 and 28 DF,  p-value: 0.0001654
## visualize model for deltametRate/avgPercLoading

par(mfrow = c(1,1))

plot(dmr_dpl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod3, newdata = ndf), ndf)

ggplot(newDF, aes(x = avgPercLoad, y = dmr_dpl)) + 
     geom_point(aes(color = order_1), shape = 17) + 
     geom_line(data = predDF, aes(x = avgPercLoad, y = preds)) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in Metabolic Rate (mL CO2 / hr) / \n Change in Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8))

ggsave("dmr_dpl.pdf", width = 5, height = 4)

(DeltaFrq^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {df2_dpl = deltaFrq2 / deltaPercLoad})
par(mfrow = c(1,1))
plot(df2_dpl ~ avgPercLoad, newDF)

mod1_f <- lm(df2_dpl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_f)
## 
## Call:
## lm(formula = df2_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.651 -21.833  -1.027  24.243  50.541 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         192.2395    27.4029   7.015 1.89e-07 ***
## avgPercLoad          -2.6608     0.4803  -5.539 8.14e-06 ***
## order_1loadedSecond -36.8089    11.4800  -3.206  0.00355 ** 
## size_pc1             -1.5020     3.0904  -0.486  0.63103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.59 on 26 degrees of freedom
## Multiple R-squared:  0.6538, Adjusted R-squared:  0.6139 
## F-statistic: 16.37 on 3 and 26 DF,  p-value: 3.534e-06
mod2_f <- update(mod1_f, .~. - size_pc1)
anova(mod1_f, mod2_f) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: df2_dpl ~ avgPercLoad + order_1 + size_pc1
## Model 2: df2_dpl ~ avgPercLoad + order_1
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     26 24327                           
## 2     27 24548 -1   -221.01 0.2362  0.631
summary(mod2_f) # not final model for paper
## 
## Call:
## lm(formula = df2_dpl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.421 -24.372   0.984  25.355  51.493 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         188.8709    26.1342   7.227 8.99e-08 ***
## avgPercLoad          -2.5918     0.4524  -5.730 4.33e-06 ***
## order_1loadedSecond -37.8192    11.1294  -3.398  0.00212 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.15 on 27 degrees of freedom
## Multiple R-squared:  0.6507, Adjusted R-squared:  0.6248 
## F-statistic: 25.15 on 2 and 27 DF,  p-value: 6.817e-07
# get p-values for avgPercLoad 
anova(mod2_f,  update(mod2_f, .~. - order_1)) # p-value for order
## Analysis of Variance Table
## 
## Model 1: df2_dpl ~ avgPercLoad + order_1
## Model 2: df2_dpl ~ avgPercLoad
##   Res.Df   RSS Df Sum of Sq      F  Pr(>F)   
## 1     27 24548                               
## 2     28 35047 -1    -10499 11.547 0.00212 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2))
plot(mod2_f) # non-linearity in residuals

plot(mod2_f, which = 4)


# update model to add a non-linear term
m22 <- lm(df2_dpl  ~ avgPercLoad_cent  + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m22)
## 
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + 
##     size_pc1, data = newDF)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -40.78 -19.93  -4.39  11.41  50.81 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          31.26529    7.86521   3.975 0.000528 ***
## avgPercLoad_cent     -2.84108    0.41899  -6.781 4.17e-07 ***
## apl_cent2             0.09699    0.03094   3.135 0.004358 ** 
## order_1loadedSecond -45.07571   10.26363  -4.392 0.000180 ***
## size_pc1             -1.13910    2.67274  -0.426 0.673614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.43 on 25 degrees of freedom
## Multiple R-squared:  0.7515, Adjusted R-squared:  0.7117 
## F-statistic:  18.9 on 4 and 25 DF,  p-value: 2.873e-07
car::vif(m22) # vif is much better with centered variables
## avgPercLoad_cent        apl_cent2          order_1         size_pc1 
##         1.135857         1.108436         1.126036         1.118753
m22a <- update(m22, .~. - size_pc1)
anova(m22, m22a) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     25 17463                           
## 2     26 17590 -1   -126.88 0.1816 0.6736
summary(m22a)
## 
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = newDF)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -42.26 -19.10  -2.71  11.28  53.26 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          31.61308    7.69867   4.106 0.000354 ***
## avgPercLoad_cent     -2.78994    0.39507  -7.062 1.69e-07 ***
## apl_cent2             0.09756    0.03042   3.207 0.003540 ** 
## order_1loadedSecond -45.88915    9.92463  -4.624 9.07e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.01 on 26 degrees of freedom
## Multiple R-squared:  0.7497, Adjusted R-squared:  0.7208 
## F-statistic: 25.96 on 3 and 26 DF,  p-value: 5.55e-08
par(mfrow = c(2,2))

plot(m22a) # residuals look better, though slight fan shape

plot(m22a, which = 4) ## row 22 looks highly influential
nd22 <- newDF[-22, ]
m22s <- lm(df2_dpl  ~ avgPercLoad_cent  + apl_cent2 + order_1, data = nd22)
summary(m22s) # no major change when we remove obs num 22
## 
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = nd22)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.038 -17.721  -6.681  11.348  54.728 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          29.57857    7.34149   4.029 0.000460 ***
## avgPercLoad_cent     -2.55219    0.39112  -6.525 7.79e-07 ***
## apl_cent2             0.12123    0.03101   3.910 0.000624 ***
## order_1loadedSecond -45.07241    9.38443  -4.803 6.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.57 on 25 degrees of freedom
## Multiple R-squared:  0.7554, Adjusted R-squared:  0.7261 
## F-statistic: 25.74 on 3 and 25 DF,  p-value: 8.179e-08
plot(m22s, which = 4)

summary(m22a) # final model for paper (though we could also report the non-centered values)
## 
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = newDF)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -42.26 -19.10  -2.71  11.28  53.26 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          31.61308    7.69867   4.106 0.000354 ***
## avgPercLoad_cent     -2.78994    0.39507  -7.062 1.69e-07 ***
## apl_cent2             0.09756    0.03042   3.207 0.003540 ** 
## order_1loadedSecond -45.88915    9.92463  -4.624 9.07e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.01 on 26 degrees of freedom
## Multiple R-squared:  0.7497, Adjusted R-squared:  0.7208 
## F-statistic: 25.96 on 3 and 26 DF,  p-value: 5.55e-08
## visualize model for deltametRate/avgPercLoading

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$apl <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)

predDF <- data.frame(preds = predict(m22a, newdata = ndf, se = TRUE), ndf)

par(mfrow = c(1,1))

plot(df2_dpl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

newDF[22, ]
##    BeeID Mstarved        S  Itspan          L order_H order_L   M2_H
## 22   E42   0.1371 6.27e-05 0.00486 0.01097648       2       1 0.2945
##      M2_L   MF_H   MF_L  MetR_H   MetR_L   freq_H   freq_L    amp_H
## 22 0.2145 0.2761 0.2002 12.2844 11.52512 175.1276 186.0056 138.0042
##       amp_L load_H  load_L   MT_H    MT_L perLoad_H perLoad_L      frce_H
## 22 117.3242 0.1482 0.07025 0.2853 0.20735  108.0963  51.23997 0.003024302
##         frce_L      arcL2_H      arcL2_L  freq2_H  freq2_L      U_H
## 22 0.002465804 0.0003931771 0.0002841704 30669.69 34598.09 6.945104
##         U_L     arcL_H     arcL_L deltaLoad2 dLoad_nonCent deltaLoad
## 22 6.271125 0.01982869 0.01685735 1.0609e-06       0.07795   0.00103
##    deltaFreq2Perc   deltaArcL2  deltaArcL deltaFrq2 deltaMetR avgPercLoad
## 22      -69.09351 0.0001090067 0.00297134 -3928.402 0.7592781    79.66813
##    deltaPercLoad  size_pc1 size_pc1_2      order_1    dmr_dpl
## 22      56.85631 0.4108091  0.1687641 loadedSecond 0.01335433
##    avgPercLoad_cent apl_cent2   df2_dpl
## 22         23.00355  529.1635 -69.09351
ggplot(newDF, aes(x = avgPercLoad, y = df2_dpl)) + 
     geom_point(aes(color = order_1),shape = 17) + 
     geom_point(aes(size = BeeID == "E42")) +  # show the influential point
     geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in wingbeat freq^2 (hz^2) / \n Change in Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8)) 
## Warning: Using size for a discrete variable is not advised.

ggplot(newDF, aes(x = avgPercLoad, y = df2_dpl)) + 
     geom_point(aes(color = order_1),shape = 17) + 
     geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in wingbeat freq^2 (hz^2) / \n Change in Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8)) 

ggsave("df2_dpl.pdf", width = 5, height = 4)

(DeltaArcL^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {da2_dpl = deltaArcL2 / deltaPercLoad})

par(mfrow = c(1,1))
plot(da2_dpl ~ avgPercLoad, ylab = c("deltaArcL^2 / avgPercLoad"), newDF)

mod1_a <- lm(da2_dpl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_a)
## 
## Call:
## lm(formula = da2_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.132e-06 -4.984e-07  6.244e-08  4.953e-07  1.054e-06 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          1.719e-06  5.312e-07   3.237  0.00329 **
## avgPercLoad          3.674e-09  9.312e-09   0.395  0.69641   
## order_1loadedSecond -2.368e-07  2.226e-07  -1.064  0.29706   
## size_pc1             4.140e-08  5.991e-08   0.691  0.49567   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.93e-07 on 26 degrees of freedom
## Multiple R-squared:  0.05092,    Adjusted R-squared:  -0.05859 
## F-statistic: 0.465 on 3 and 26 DF,  p-value: 0.7092
mod2_a <- update(mod1_a, .~. - avgPercLoad)  # I also checked for a squared term in this model (code not shown)
anova(mod1_a, mod2_a) # remove avgPercLoad
## Analysis of Variance Table
## 
## Model 1: da2_dpl ~ avgPercLoad + order_1 + size_pc1
## Model 2: da2_dpl ~ order_1 + size_pc1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     26 9.1429e-12                             
## 2     27 9.1976e-12 -1 -5.4736e-14 0.1557 0.6964
summary(mod2_a)
## 
## Call:
## lm(formula = da2_dpl ~ order_1 + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.103e-06 -5.079e-07  7.199e-08  4.954e-07  9.781e-07 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.919e-06  1.568e-07  12.242 1.57e-12 ***
## order_1loadedSecond -2.214e-07  2.156e-07  -1.027    0.314    
## size_pc1             3.442e-08  5.634e-08   0.611    0.546    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.837e-07 on 27 degrees of freedom
## Multiple R-squared:  0.04524,    Adjusted R-squared:  -0.02549 
## F-statistic: 0.6396 on 2 and 27 DF,  p-value: 0.5353
mod2_b <- update(mod2_a, .~. - size_pc1)
anova(mod2_a, mod2_b)  # remove size
## Analysis of Variance Table
## 
## Model 1: da2_dpl ~ order_1 + size_pc1
## Model 2: da2_dpl ~ order_1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     27 9.1976e-12                             
## 2     28 9.3248e-12 -1 -1.2716e-13 0.3733 0.5463
mod2c <- update(mod2_b, .~. - order_1)
anova(mod2c, mod2_b) # remove order
## Analysis of Variance Table
## 
## Model 1: da2_dpl ~ 1
## Model 2: da2_dpl ~ order_1
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1     29 9.6334e-12                            
## 2     28 9.3248e-12  1 3.0861e-13 0.9267  0.344
summary(mod2c) # final mod for paper
## 
## Call:
## lm(formula = da2_dpl ~ 1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.206e-06 -4.613e-07 -1.433e-08  4.925e-07  9.636e-07 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.801e-06  1.052e-07   17.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.764e-07 on 29 degrees of freedom
par(mfrow = c(2,2))
plot(mod2_a)

plot(mod2_a, which = 4)

# visualize model

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)

predDF <- data.frame(preds = predict(mod2c, newdata = ndf, se = TRUE), ndf)

par(mfrow = c(1,1))

plot(da2_dpl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ggplot(newDF, aes(x = avgPercLoad, y = da2_dpl)) + 
     geom_point(aes(color = order_1), shape = 17) + 
     geom_line(data = predDF, aes(x = avgPercLoad, y = preds.fit)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in arc length^2 (radians^2) / \n Change in Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8)) 

ggsave("da2_dpl.pdf", width = 5, height = 4)

Session Info

sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plyr_1.8.4         viridis_0.3.4      visreg_2.3-0      
##  [4] effects_3.1-2      sjPlot_2.2.1       influence.ME_0.9-8
##  [7] reshape2_1.4.2     MASS_7.3-45        gsheet_0.4.2      
## [10] lme4_1.1-12        Matrix_1.2-8       car_2.1-4         
## [13] ggplot2_2.2.1     
## 
## loaded via a namespace (and not attached):
##  [1] tidyr_0.6.1        splines_3.3.2      merTools_0.3.0    
##  [4] modelr_0.1.0       shiny_1.0.0        assertthat_0.1    
##  [7] stats4_3.3.2       coin_1.1-3         yaml_2.1.14       
## [10] backports_1.0.5    lattice_0.20-34    quantreg_5.29     
## [13] digest_0.6.12      minqa_1.2.4        colorspace_1.3-2  
## [16] sandwich_2.3-4     htmltools_0.3.5    httpuv_1.3.3      
## [19] psych_1.6.12       broom_0.4.1        SparseM_1.74      
## [22] haven_1.0.0        purrr_0.2.2        xtable_1.8-2      
## [25] mvtnorm_1.0-5      scales_0.4.1       stringdist_0.9.4.4
## [28] MatrixModels_0.4-1 arm_1.9-3          tibble_1.2        
## [31] mgcv_1.8-16        DT_0.2             TH.data_1.0-8     
## [34] nnet_7.3-12        lazyeval_0.2.0     pbkrtest_0.4-6    
## [37] mnormt_1.5-5       survival_2.40-1    magrittr_1.5      
## [40] mime_0.5           evaluate_0.10      nlme_3.1-130      
## [43] foreign_0.8-67     tools_3.3.2        multcomp_1.4-6    
## [46] stringr_1.1.0      munsell_0.4.3      blme_1.0-4        
## [49] grid_3.3.2         nloptr_1.0.4       htmlwidgets_0.8   
## [52] labeling_0.3       rmarkdown_1.3      gtable_0.2.0      
## [55] codetools_0.2-15   curl_2.3           abind_1.4-5       
## [58] sjstats_0.7.1      DBI_0.5-1          sjmisc_2.2.1      
## [61] R6_2.2.0           gridExtra_2.2.1    zoo_1.7-14        
## [64] knitr_1.15.1       dplyr_0.5.0        rprojroot_1.2     
## [67] readr_1.0.0        modeltools_0.2-21  stringi_1.1.2     
## [70] parallel_3.3.2     Rcpp_0.12.9        coda_0.19-1       
## [73] lmtest_0.9-34
Sys.time()
## [1] "2017-03-16 17:20:54 EDT"